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We calculate the equation of state (EOS) of dense matter, using a relativistic mean field (RMF) 
model with a density dependent coupling that is a slightly modified form of the original NL3 interac- 
tion. For nonuniform nuclear matter we approximate the unit lattice as a spherical Wigner-Seitz cell, 
wherein the meson mean fields and nucleon Dirac wave functions are solved fully self-consistently. 
We also calculate uniform nuclear matter for a wide range of temperatures, densities, and proton 
fractions, and match them to non-uniform matter as the density decreases. The calculations took 
over 6,000 CPU days in Indiana University's supercomputer clusters. We tabulate the resulting 
EOS at over 107,000 grid points in the proton fraction range Yp = to 0.56. For the temperature 
range T — 0.16 to 15.8 MeV we cover the density range ub = 10^* to 1.6 fm^''; and for the higher 
temperature range T = 15.8 to 80 MeV we cover the larger density range ub ~ 10^** to 1.6 fm~'^. 
In the future we plan to study low density, low temperature (T<15.8 MeV), nuclear matter using a 
Virial expansion, and we will match the low density and high density results to generate a complete 
EOS table for use in astrophysical simulations of supernova and neutron star mergers. 
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I. INTRODUCTION 



The equation of state (EOS) for hot, dense matter in 
massive stars relates energy and pressure to temperature, 
density, and composition. It has been a long-standing 
problem to understand the EOS at both subnuclear and 
supranuclear density, to which great efforts have been 
devoted, from laboratory heavy ion colhsion experiments 
[H , computer simulations of supernova , and theo- 
retical many-body calculations [J| . The EOS of hot dense 
matter in supernovae (SN) and neutron star (NS) merg- 
ers encompass multi-scale physics. Temperature can vary 
from to as high as 100 MeV, exciting nuclei, nucleon 
and possibly pion and other degrees of freedom. The den- 
sity can vary from k, 10^ to 10^^ g- cm~'^, where matter 
can be in gas, liquid or solid phases. The proton fraction 
can vary from to 0.6, from extremely neutron rich mat- 
ter to proton rich matter. These very large parameter 
ranges make construction of a full EOS table difficult. It 
is necessary to employ different approximations for dif- 
ferent parameter ranges. As a result, there exist only two 
realistic EOS tables that are in widespread use for astro- 
physical simulations, the Lattimer-Swesty (L-S) equation 
of state @, that uses a compressible liquid drop model 
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with a Skyrme force, and the H. Shen, Toki, Oyamatsu 
and Sumiyoshi (S-S) equation of state [1, 0, that uses 
the Thomas Fermi and variational approximations with 
a relativistic mean field (RMF) model. We plan to gen- 
erate a complete equation of state, employing relativistic 
mean field calculations for matter at intermediate and 
high density as described in this paper. In the future 
we plan to use the Virial expansion of a nonideal gas to 
describe matter at low density. These two parts will be 
matched together and we will generate a thermodynam- 
ically consistent EOS over the full range of parameters. 
Finally we will generate additional EOSs from RMF mod- 
els with different high density symmetry energies. This 
will allow one to correlate features of astrophysical simu- 
lations with properties of the symmetry energy assumed 
for the EOS. 

There are still large uncertainties in the EOS at 
supranuclear densities. The density dependence of the 
symmetry energy dS/driB is poorly known and strongly 
influences the stiffness of the EOS. It can be constrained 
from measurements of NS radii and masses Q , precision 
determination of the neutron rms radius in ^°^Pb 0], and 
also heavy ion collision experiments [l| . A stiff EOS (high 
pressure) at high density gives larger NS radii, while a 
stiff EOS at normal and low density favors a larger neu- 
tron radius in ^ospj^ rpj^g elliptic and transverse 
flow observables in heavy ion collisions are sensitive to 
the isospin dependence of mean fields and to the EOS 
at densities up to a few times nuclear saturation den- 
sity. Many nuclear many-body models fall into two cat- 
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egories, the non-relativistic Skyrme models (See for ex- 
ample, Ref. [ll[ for a review) and relativistic mean field 
models ■ The parameters in these models are usually 
fitted to nuclear properties at normal nuclear densities, 
afterwards they are extrapolated to study supranuclear 
matter. The L-S EOS uses a Skyrme model featuring 
a relatively soft EOS and the S-S EOS uses the RMF 
interaction TMl that features a stiffer EOS. Since the 
symmetry energy is not well constrained, it is important 
to explore the effects of different symmetry energies on 
the EOS and SN simulations. 

In this paper, we use a RMF model for non-uniform 
matter at intermediate density and uniform matter at 
high density. Low density pure neutron matter is anal- 
ogous to a unitary gas [l3|, where the neutron- neutron 
scattering length is much larger than both the effective 
range and the average inter-particle spacing. To bet- 
ter describe neutron-rich matter at low density, we use 
a density dependent scalar meson-nucleon coupling. At 
high density, the model reduces to the normal RMF pa- 
rameter set NL3. The unit lattice of non-uniform nu- 
clear matter is conveniently approximated by a spherical 
Wigner-Seitz (W-S) cell. The meson mean fields and nu- 
cleon Dirac wave functions inside the Wigner-Seitz cell 
are solved fully sclf-consistcntly. This is unlike the S-S 
EOS that used Thomas-Fermi and variational approxi- 
mations and the L-S EOS that used a simple liquid drop 
model. The size of the W-S cell is found by minimization 
of the free energy per nucleon. The W-S approximation 
provides a framework to incorporate the best known mi- 
croscopic nuclear physics • The nuclear shell structure 
effects are included automatically and it is already pos- 
sible for some effects of complex nuclear pasta states to 
be included in spherical calculations in the form of shell 
states [15|. Full three-dimensional W-S calculations in 
principle could incorporate various pasta shapes [16, 17], 
which would make the transition to uniform matter more 
smooth. However this will demand much larger compu- 
tational resources. In this work we use the spherical W-S 
approximation. 

Our relativistic mean field calculations can accurately 
describe the radial shape of large neutron rich nuclei in- 
cluding the expected neutron rich skin. In contrast the 
original L-S EOS is based on a very simple liquid drop 
model of nuclear structure that may incorrectly describe 
the neutron skin. Alternatively the S-S EOS is based on a 
Thomas Fermi approximation that neglects shell effects. 
These are included in our Hartree calculations. Further- 
more, the variational forms for the densities assumed 
by S-S may be a poor approximation for large proton 
numbers where the Coulomb repulsion is large. Instead 
our exact solutions of the radial mean field equations 
allow richer density distributions including shell states 
with central depressions [l5|. These differences in densi- 
ties may be important for neutrino interactions in Super- 
novae. Finally our calculations correctly reproduce the 
Unitary gas limit for a low density neutron gas, see be- 
low, while both the L-S and S-S EOS reduce incorrectly 



to the energy of free neutrons. 

One can demand that any EOS be consistent with, pos- 
sibly model dependent interpretations of, observations of 
neutron stars. For example, Klahn et. al. propose a se- 
ries of tests that an EOS should satisfy to be consistent 
with observations [l^. They demand that any reliable 
nuclear EOS be able to reproduce the recently reported 
high pulsar mass of 2.1± 0.2 Mq for PSR J0751-H1807 
18 1 . However, this observation may have been retracted 
20|. Furthermore, Klahn etal. require the EOS to repro- 



duce a large binding energy for Pulsar B in J0737-3039. 
However, this conclusion could be sensitive to assump- 
tions about the system such as the amount of mass loss. 
Klahn et. al. go on to demand that the EOS not allow 
direct URCA cooling of neutron stars of mass 1 to 1.5 
Mq. We consider a more conservative approach. While 
many stars cool slowly, observations do suggest that at 
least some stars have enhanced cooling. Unfortunately 
observations do not directly constrain the mass that may 
separate enhanced from normal cooling. Indeed, there 
is little direct observational evidence that more massive 
stars cool more quickly, although this is a theoretical prej- 
udice. 

One can also use laboratory data to constrain the EOS. 
The neutron skin thickness of a heavy nucleus constrains 
the density dependence of the symmetry energy. Further- 
more, there are many measurements of the skin thickness 
with a variety of strongly interacting probes. However, 
there may be important model dependence from strong 
interaction uncertainties. For example (^He,t) measure- 
ments of spin dipole strength have been used to extract 
neutron skin thicknesses in Sn isotopes [2l|. For these 
measurements, the spin dipole strength was assumed to 
be proportional to the measured cross section, and the 
proportionality constant was arbitrarily fixed in order to 
reproduce the skin thickness of -'^^"Sn as predicted by an 
old Hartree- Fock calculation Presumably, if a dif- 

ferent skin thickness in "'^^"Sn is fit, one would also get 
different answers for the skin thickness in other isotopes. 

This situation may soon change. The Lead Radius Ex- 
periment (PREX) at Jefferson Laboratory is using par- 
ity violating electron scattering to measure the neutron 
skin thickness in ^°^Pb Q. Parity violation is a sensitive 
probe of neutrons because the weak charge of a neutron 
is much larger than that of a proton. Furthermore, this 
electro-weak reaction may have much smaller strong in- 
teraction uncertainties. Data taking for PREX should be 
completed by June 2010. 

Instead of trying to determine, ahead of time, the 
best EOS to satisfy existing observational constraints, we 
adopt what we hope will be a more robust approach. We 
are calculating a number of EOSs based on different ef- 
fective interactions. In this paper we present first results 
for the NL3 interaction with a symmetry energy that is 
large at high densities. In later work we will present 
EOSs with softer high density symmetry energies. These 
different EOS will allow one to correlate features of as- 
trophysical simulations with properties of the EOS. Then 



3 



one can draw conclusions based on combined information 
from laboratory experiments and astronomical observa- 
tions. 

In this paper we focus on nucleon degrees of freedom. 
Hyperons could play a role at high densities, see for ex- 
ample (23I I. However, the contribution of hyperons could 
depend on uncertain hyperon interactions. In addition, 
there could be pion or kaon condensates or a variety of 
quark matter phases. See for example the review by Page 
and Reddy [3| . Chiral symmetry restoration and the soft- 
ening of pionic or kaonic modes could be important. Fi- 
nally, thermal pions and pion interactions should be very 
important at high temperatures. All of these effects may 
increase the uncertainties in the EOS. 

We tabulate the equation of state for intermediate and 
high density nuclear matter over the range of tempera- 
tures T, densities ub, and proton fractions Yp given in 
Table U and described in Sec. IIVI We calculate the free 
energy of nonuniform matter for 17021 points, and the 
free energy of uniform matter for 90478 points in T, n^, 
and Yp space. This took 6000 CPU days on Indiana 
University's supercomputer clusters. 

The paper is organized as follows: in Section |ll] the 
density dependent RMF model is explained in detail. In 
Section Hm we describe the RMF parameters that we use 
including a density dependent coupling. We describe the 
computational methodology for our large computer runs 
in Section ITVl Section fVl shows results for RMF calcula- 
tions, including the free energy and the nucleon density 
distributions in the non-uniform W-S cells. Finally, Sec- 
tion I VII presents a summary of our results and gives an 
outlook for future work. 



TABLE I: Range of temperatures T, densities ns, and proton 
fractions Yp in the EOS table. 



Parameter 


Low T 


High T 


Total # 


logio(T) [MeV] 


-0.8 to 1.2 


1.2 to 1.9 


32 


logio(ns) [fm^^] 


-4.0 to 0.2 


-8.0 to 0.2 


43, 83 


Yp 


0,0.05 to 0.56 


0,0.05 to 0.56 


53 



II. FORMALISM 

We now describe the mean field formalism that we use 
for non-uniform matter in Section III Al and for uniform 
matter in Section fllBI 



A. Non-uniform nuclear matter in Wigner-Seitz 
approximation 

The formalism for relativistic mean field theory has 
been reviewed in previous papers, see eg [l^]- To better 
describe neutron rich matter at low density we introduce 
a density dependent coupling between the scalar meson 
and the nucleon as described in Section illll We note that 



many previous studies of density dependent RMF models 
mainly focused on better descriptions of nuclear matter 
at supranuclear density (see for example, [H, H^). In 
this section we focus on low density neutron rich matter. 

The basic ansatz of the RMF theory is a Lagrangian 
density where nucleons interact via the exchange of 



sigma- (cr), omega- (i:^^) 
photons {Afj). 



and rho- (p^) mesons, and also 



1 + ^3 
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We note that T^r — Po- (n) {n = ^f]^ and is nucleonic 
current) is the density dependent coupling between the 
sigma meson and the nucleon. Here the field tensors of 
the vector mesons and the electromagnetic field take the 
following forms: 



9pP X P 



(2) 



In charge neutral nuclear matter composed of neu- 
trons, n, protons, p, and electrons, e, there are equal 
numbers of electrons and protons. Electrons can be 
treated as a uniform Fermi gas at high densities [33 |. 
They contribute to the Coulomb energy of the npe mat- 
ter and serve as one source of the Coulomb potential. 

The variational principle leads to the following equa- 
tions of motion 

[a-p + V{r) + I3{m + S'(r))]V'i = SiV'i (3) 

for the nucleon spinors, with vector and scalar potentials 



V{r) = + .g,f • + A, 

^(r) = P.a, 



where 



n dn 



(4) 



(5) 



is the rearrangement term due to the density dependent 
coupling between the sigma meson and the nucleon, and 
Ps is the scalar density of nucleons to be defined below. 

The equations of motion for the mesons and photons 
are, 



{ml - V^) cr = -T^ps - 92(7^ - 53C^^ 
{ml - V2)w^= gjf^ -c^uJ^^{uJ-uJ,), 
{ml - V2)pf = g,Jf , 

-^^A^^= e{j^-j^:), 



(6) 
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where the electrons are included as a source of Coulomb 
potential. The nucleon spinors provide the relevant 
source terms: 

(7) 

At finite temperature, Fermi-Dirac statistics imply the 
occupations, n^, of protons and neutrons are: 



where /i is the chemical potential for neutron (proton). 
In the calculation, we include all levels with gi-ni > 10"^, 
where gi is the degeneracy of the level. 

Since the systems under consideration have tempera- 
tures of, at most, tens of MeV, we neglect the contribu- 
tion of negative energy states, i.e., the so-called no sea 
approximation. In a spherical nucleus, there are no cur- 
rents in the nucleus and the spatial vector components 



of w^, Pfj, and A^^ vanish. One is left with the timelike 
components, uiq, po and Aq. Charge conservation guar- 
antees that only the 3-component of the isovector po.s 
field survives. The above non-linear equations are solved 
by iteration within the context of the mean field approx- 
imation whereby the meson field operators are replaced 
by their expectation values. 

The spherical Wigner-Seitz approximation is used to 
describe non-uniform matter. One W-S cell has one nu- 
cleus. In this approximation it is important to include 
lattice Coulomb corrections between neighboring W-S 
cells. The detailed treatments we use have been discussed 
in a previous paper [l5l | and we will not repeat them here. 

We solve for the meson mean fields and the nucleon 
Dirac wave functions self-consistently inside a W-S cell of 
radius Rc, for a given baryon density ub, proton fraction 
Yp and temperature T. In our RMF model nucleons 
(proton and neutron) are the only baryons. The nucleon 
number inside the W-S cell is ^ = AttR^ub and the 
proton number is Z ^ YpA. The internal energy of a W- 
S cell, including the approximate lattice Coulomb energy 
correction, is. 



Eb — Enucleon + + Ep + + Ecoul — mA 



/dV If 11 

d^Tja{r)-^Ps{r)a{r) - ^ J d^f{^<T(TPs{r) + -g2(T^ + ^Ss^"}, 

\J '^^'^9pPo,3jo,3{^) ~ \J d^'^iduji^o joir) - ^caWo}- ^ J {pp + Pe)Ao{r)d^r + dw - mA, (9) 



where dw = 0.0065620Z2/a is the approximate Coulomb 
correction for a bcc lattice(2^], and = Vw-s is the 
volume of W-S cell. 

The nucleon contribution to the entropy is given by 
the usual formula, 

Sb = -kB^gi[niln{ni) + (1 - ni)ln(l - n^)] , (10) 

i 

where Ui is given in Eq. (g]). With Eqs. dH) and ([TU]). 
it is easy to obtain the nucleon contribution to the free 
energy per nucleon F, 

F = Fb/A = {Eb - TSb)/A. (11) 

B. Uniform nuclear matter 

To make the paper self-contained, we give the formulas 
for uniform matter in RMF model. As we show below. 



at high temperatures or high densities the matter is uni- 
form. We include anti-nucleon terms which make a small 
contribution at very high temperatures. 

The energy density of uniform nuclear matter is, 

i=N,P 

+ ^52^^ + \g3^^ + jcgW^, (12) 

where 

4» = / d^kE*{k)[nk{T)+n^{T)], (13) 

with effective mass m* — m + T^a, and E*{k) — 
Vfc^ + ra*"^ . The occupation probabilities for particles 
nk{T) and antiparticles nk{T) are, 
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nk{T) 
nk{T) 



1 

eK-p{E* (k) + g^ojo + gpT^PQ,3 + ^Ps<J - p)/T + 1 ' 

1 

exp(£'*(fc) - guiUjQ - gpT^po,3 - ^Pscr + p)/T + 1 



(14) 
(15) 



The pressure of uniform nuclear matter is, 

i=N,P 

(16) 

I 



where 



kin 



d-^k- 



[nk{T)+nk{T)]. (17) 



3(27r)3 

The entropy density of uniform nuclear matter is, 



2kB 
(27r)3 



Sk[nk{T)\nnk(T) + {I - nk{T))\n{l - nk{T)) + nk{T)\nnk{T) + (1 - nfe(T))ln(l - nfe(T))].(18) 

I 



Using Eqs. ([T2|) and (fTSl) one can obtain the free energy 
density per nucleon for uniform matter, 



F 



(e - Ts)/nB- 



(19) 



III. 



PARAMETER SET WITH DENSITY 
DEPENDENT COUPLING 



In this work we use the NL3 effective interaction [23] 
that has been successful in reproducing ground state 
properties of stable nuclei and the saturation properties 
of symmetric nuclear matter. The values of parameters 
in the NL3 effective interaction are listed in Table HIl 



TABLE II: NL3 effective interaction. The nucleon masses are 
M — 939 MeV for both protons and neutrons and C3 = in 
Eq.p. 



9p 



92 
(fm-i) 



53 rn^ m,ui mp 
(MeV) (MeV) (MeV) 



10.217 12.868 4.474 -10.431 -28.885 508.194 782.5 763 

As is well known, the mean field approach for pure 
neutron matter is problematic at low densities because 
long range correlations are important. Neutron mat- 
ter at low density is very close to a unitary gas [T3| . 
since the scattering length is much larger than the inter- 
particle spacing, which is also larger than the effective 
range of nuclear interaction. To describe neutron mat- 
ter phenomenologically in the RMF framework, without 
losing its success for the properties of nuclear matter, 
we introduce a density-dependent scalar meson-nucleon 



coupling. 



l + Q 



n > no 
n < no- 



(20) 



The two free parameters hq and a are determined by 
matching the energy of neutron matter to that of a uni- 
tary gas at zero temperature Ejj 29], 



El 



3 k 



3fc 



^.--JL^O,U---^, (21) 



5 2m 



5 2m 



where fc^^^ is the neutron Fermi momentum. The best 
fitted values are no = 5 x 10"'^ fm~'^ and a = 1.2. 

In Fig.m the energy of pure neutron matter at T = is 
shown for the original NL3 set, the modified NL3 set with 
a density dependent cr-N coupling as in Eq. (PH)) . and the 
unitary gas calculated by Eq. (PT|) . The unitary gas gives 
lower energy than the original NL3 result by about 0.2 
MeV per particle, due to the strong S wave attractive 
interactions. This energy difference is very relevant for 
matching a Virial expansion to the mean field calcula- 
tions, since the Virial expansion includes the lon g ra nge 
two-body neutron-neutron attractive interaction [23, |30] 
while the normal mean field calculation does not. In the 
density range shown in the figure, NL3 also gives a lower 
energy than the TMl or FSUGold [U RMF parameter 
sets. However, the density dependent NL3 can fit the 
unitary gas result by tuning the coupling strength in the 
attractive scalar meson channel. Therefore, the density 
dependent NL3 set describes successfully the properties 
of both neutron rich matter and nuclear matter. In the 
following when we refer to NL3 set, we mean the density 
dependent NL3 set unless otherwise specified. 
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lO"*^ 10"^ 10* 10"^ 10"^ 

n [fm" ] 



FIG. 1: (Color on line) Energy of pure neutron matter at T 
= 0. The red curve is from the original NL3 set. The black 
curve is for NL3 with a density dependent cr-N coupling To- 
as in Eq. (|20p . The blue dashed line is the energy of a unitary 
gas, see Eq. (pTjl . 



IV. COMPUTATIONAL METHODOLOGY 

In this section we describe our strategy for evaluating 
the equation of state. We calculate the equation of state 
for the following partitioning of T, n^, and Yp parameter 
space, see Tabic HI 

• We use a step of 0.2 in logio(r/ [MeV]) for 
logio{T/ [MeV]) from -0.8 to 0, a step of 0.1 for 
logio(T/[MeV]) from to 1.1 and a step of 0.05 
for login (T/ [MeV]) above 1.1. We have a total of 
32 points for T from 0.16 to 80 MeV. 

• For temperatures T below 15.8 MeV (where mat- 
ter can be nonuniform) we use a step of 0.1 in 
logio(nB/ [fm-3]) for logiglns/ [fm"^]) from -4.0 
to 0.2. We have a total of 43 points for ns from 
10-4 to 1.6 fm-3. 

• For temperatures T above 15.8 MeV (where 
matter is uniform) we use a step of 0.1 in 
logio(riB/ [fm-3]) for logioCns/ [fm'^j ) from -8.0 
to 0.2. We have a total of 83 points for ub from 
10-^ to 1.6 fm-3. 

• We use a step of 0.01 in proton fraction Yp for Yp 
from 0.05 to 0.56. We aslo include Yp = 0.0. This 
gives a total of 53 points for Yp from 0.0 to 0.56. 

This partitioning gives a total of 40,248 points in the 
nonuniform Hartree region. However for matter at higher 
temperatures, but still T < 15.8 MeV, and/or lower 
proton fractions, Hartree results give higher free ener- 
gies than corresponding results for uniform matter. By 
roughly estimating the phase boundary, and keeping 
enough points to cross the transition density, see Section 



IVl , we calculate the free energy for a reduced number of 
points in the nonuniform Hartree region, see Sec. Ill A[ 
which includes a total of 17,021 points. We also calcu- 
late free energies for uniform nuclear matter at a total of 
90,478 points, see Sec. IiTbI 

The most time is spent evaluating (temperature T, 
proton fraction Yp, density ns) points in the non- 
uniform Hartree mean field region. For each point we 
need to minimize the free energy of the W-S cell with 
respect to the cell radius which typically requires evalu- 
ation at 40 to 100 cell radii. This minimization can be 
complicated by the existance of local minima. For each 
cell size, we need to solve the mean fields self consistently. 
We have already developed the code for this minimiza- 
tion (see Ref. jT^) which is slightly modified in this 
work to accomodate the density dependent coupling in 
RMF. 

The mean fields provide potentials for the individual 
nucleons in the W-S cell which obey the Dirac equation. 
The Dirac equation is solved by a fourth order Runge- 
Kutta method with shooting techniques. For nuclear 
matter at finite temperature, there could be hundreds of 
nucleons that populate thousands of levels according to 
Fermi-Dirac statistics. For each level, we need to solve 
the Dirac equation. The potentials for the nucleons in 
the Dirac Eq. ^ are various meson mean fields which 
obey the extended Klein-Gordon (K-G) equation. The 
source terms for the K-G equations are provided by var- 
ious nucleon density terms in Eq. ([7]). Given the nucleon 
density terms, the K-G equations are solved by a Green's 
function method, which updates the meson mean fields. 
The updated mean fields can now be used to solve the 
nuclear levels and nucleon densities again. This process 
is repeated until full self-consistency is achieved in both 
the mean fields and the nuclear levels. 

Computationally, the problem to be solved is embar- 
rassingly parallel because each point of density, temper- 
ature and proton fraction is independent of the others. 
A total number of ^ 17,000 independent tasks must be 
run, where each task calculates the required quantities 
at a single point in the phase space. Unfortunately, the 
run time on an individual task varies from a few minutes 
to more than 24 hours, depending on the number of iter- 
ated cell radii, and the number of nucleon energy levels 
included. 

Each point in the phase space was mapped to a 
unique integer that we refer to as the job index. A 
file, runlist, was prepared with a list of job indi- 
cies for the whole phase space, and a single character 
(A=available, R=running, r=Re-running, C=complete, 
T=time-limited and F=failed) that gives the status of 
calculations for that job index. An MPI parallel wrapper 
code manages the running of the many requested tasks. 
Typically, one parallel job requests a set of compute cores 
(usually 256). Each MPI rank, using a single CPU core, 
is assigned one job index corresponding to one point in 
the phase space and it evaluates the required quantities. 

Initially, rank zero of the MPI job 
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• locks the job listing file runlist 

• reads runlist until a list of available tasks is filled 

• closes runlist and releases the lock 

• passes a job index to each MPI rank and begins the 
calculation for that job index. 

When the calculation completes (or time-limits or fails) 
for a given MPI rank, the status character for the job in- 
dex in runlist is modified appropriately. The now avail- 
able MPI rank will search runlist for next available task 
and the calculation restarts for the new job index. Since 
completion occurs asynchronously file locking is not used 
for this part of the process. 

A simple batch job runs through the points in phase 
space. A wall clock limit (48 hours) larger than the aver- 
age run time is used. Each rank of the MPI job can run 
a series of points via above procedure, efficiently using 
each available core for the requested wall clock period. 
One job per core is running when the wall clock limit is 
reached. These jobs are identified by being left in the " R" 
state after the batch job completes. Using this scheme 
we have achieved 85% efficiency in CPU usage. Specifi- 
cally, 85% of all jobs ended in the "C" state rather than 
the " T" or " R" state. After the runlist has been searched 
once, the remaining jobs have "R" or "T" state. Then 
these remaining jobs are resubmitted via the MPI wrap- 
per code, requesting longer time limit (typically 7 days) 
but fewer CPU cores. This procedure allows us to calcu- 
late > 99 % of the points in the runlist file. 



V. RESULTS 

In this section we discuss our results for various re- 
gions of parameter space. First, the uniform matter EOS 
at zero temperature is presented. Second, we discuss 
the free energy per nucleon for mean field calculations of 
nonuniform matter. Finally, we show the density distri- 
butions of neutrons and protons inside W-S cells. 
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FIG. 2: (Color on line) Energy of uniform matter at zero 
temperature with different proton fractions, Yp = 0, 0.1, 0.2, 
0.3, 0.4 and 0.5. 



B. Free energy and phase boundaries 

In Fig. [31 the free energy per nucleon F j A is shown 
as a function of density ub at T = 1, 3.16, 6.31 and 10 
MeV. At intermediate densities, F j A is calculated from 
Eq. (Ilip for W-S cells using Hartree mean field calcula- 
tions. At high densities, F j A is calculated from Eq. 
for uniform matter. The transition (as the density grows) 
is found at the density where uniform matter gives a lower 
free energy. In each panel, the (red) solid curves give the 
transition densities to uniform matter. The transition 
densities increase as proton fraction grows. Non-uniform 
matter can exist until higher densities in more symmet- 
ric nuclear matter. At density around 0.16 fm~^, there 
is always a minimum in the free energy per nucleon, as 
long as the proton fraction is not too small and tempera- 
ture not too high. This is the manifestation of saturation 
density in nuclear matter. 



A. Uniform matter at zero temperature 

Figure [5] shows the equation of state of uniform matter 
at zero temperature with different proton fractions Yp 
= 0, 0.1, 0.2, 0.3, 0.4, and 0.5. The solid curves are for 
the NL3 parameter set, which is used in our RMF cal- 
culations. The dashed curves show results for the TMl 
interaction, which is used in the equation of state ob- 
tained by H. Shen et al Jj]. The two sets agree to a great 
extent for densities below 0.2~0.25 fm^"^, depending the 
value of Yp. Above these densities, NL3 gives a much 
stiffer equation of state for uniform matter. This serves 
as one motivation for our choice of the NL3 parameters: 
to explore the equation of state with a stiffer symmetry 
energy. 



C. Density distributions inside Wigner-Seitz cells 

The Hartree mean field calculation provides detailed 
wavefunctions for nucleons in the non-uniform phase. In 
our W-S approximation at intermediate densities, we find 
a "spherical pasta" phase where the proton density dis- 
tribution forms a shell state with a reduced density in the 
center. This reduces the large coulomb repulsion between 
protons and was discussed in [l^ . In this section we dis- 
cuss density distributions inside the W-S cells, both for 
normal nuclei and for these shell states. 

In Fig. m neutron and proton distributions, inside the 
W-S cell, are shown for four different baryon densities, 
with r = 1 MeV and Yp — 0.1. At very low density ub 
= 0.002 fm~'^, the Hartree calculation has a minimum 
for Z = 39 protons and A — 390 nucleons. Most of 
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the neutrons are located within 10 fm of the cell center, 
although the W-S cell radius is around 31 fm. A small 
fraction of the neutrons extend to the edge of cell since 
this is an extremely neutron rich system. As the density 
rises to 0.02 fm^^^ ^j^g ^^^^ j^^g Z = 42, A = 420, 
and the neutron density at large r becomes much greater. 
The W-S cell radius drops to 17.5 fm because the lattice 
becomes more closely packed as the density increases. 
At a density of 0.05 fm'^, the W-S cell has Z = 413 
and A = 4130 and forms a shell state with both inside 
and outside surfaces. This has been discussed in our 
early paper [l^. As a result, the W-S cell radius become 
larger. At the higher density of 0.063 fm~'^, the system 
becomes uniform. 

In Fig. O the distribution of neutrons and protons are 
shown for — 0.020 fm^'^ and proton fraction Yp =0.3. 
At low temperature T = 1 MeV, where Z = 85, A = 282, 
the density distributions are similar to those for normal 
isolated nuclei. As the temperature rises to 3.16 and 
6.31 MeV, the size of W-S cell remains nearly fixed, but 
the neutron density increases with temperature at large 
radius. This is due to excitation of states with high an- 
gular momentum and/or large main quantum number as 
the temperature rises. When the temperature rises to 10 
MeV, the proton density also rises at large r, accompa- 
nied by an increase of the W-S cell size, with Z = 123, 
A = 410. At even higher temperature, the nucleus melts 
and uniform matter appears. 

Similar to Fig. [5] but at higher ub = 0.050 fm^'^ and 
proton fraction Yp = 0.45, the density distribution of 
neutrons and protons are shown for four different tem- 
peratures in Fig. [6] Here a shell state exists up to high 
temperatures. At low temperatures T = 1, 3.16 MeV, 
Z = 1315, A — 2922 and the shell state has inside and 
outside voids. As the temperature rises, nucleons pop- 
ulate both the inside and outside voids, due to thermal 
excitations. Finally, the size of the shell state shrinks at 
high temperature so that Z = 648, A = 1440 at T =10 
MeV. 



VI. SUMMARY AND OUTLOOK 



In this paper we present large scale relativistic mean 
field calculations for nuclear matter at intermediate and 
high densities. We use a density dependent modification 
of the NL3 interaction in a spherical Wigner-Seitz ap- 
proximation. Nuclear shell effects are included. We cal- 
culate the free energy, and tabulate the resulting equation 
of state at over 107,000 grid points in the proton fraction 
range Yp = to 0.56. For low temperatures T — 0.16 to 
15.8 MeV we calculate for the density range — 10"'' 
to 1.6 fm-3. For high temperatures T = 15.8 to 80 MeV, 
where the matter is uniform, we calculate for the larger 
density range ub = 10~^ to 1.6 fm~'^. These calculations 
took over 6000 CPU days. 

We solve for the nucleon Dirac wave functions and me- 
son mean fields self-consistently. This allows us to study 
how the distribution of neutrons and protons inside a 
Wigner Seitz cell evolve with density and temperature. 
We find a large variety of possible sizes and shapes for 
these distributions. 

This paper provides part of our results for an equation 
of state which will cover a broad range of temperatures, 
densities, and proton fractions. In the future, we plan 
to study low density nuclear matter using a Virial ex- 
pansion for a non-ideal gas consisting of nucleons and 
thousands of species of nuclei. Then, we will generate a 
complete thermodynamically consistent equation of state 
by matching the low density and higher density results. 
This equation of state avoids the Thomas Fermi and vari- 
ational approximations of the H. Shen et al. equation of 
state and is exact in the low density limit. It can be used 
in supernova and neutron star merger simulations. Fi- 
nally, in future work we will generate equations of state 
using other modern relativistic mean field interactions 
such as FSUGold [U. 
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FIG. 3: (Color on line) Free energy per nucleon of nuclear matter at different temperature and proton fractions. 
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FIG. 4: (Color on line) Density distribution of neutrons, protons inside W-S cell for four different baryon densities at T = 1 
MeV and proton fraction Yp = 0.1. 
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FIG. 5: (Color on line) Density distribution of neutrons, protons inside W-S cell for four different temperatures with ub = 
0.020 fm~^ and proton fraction Yp = 0.3. 
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FIG. 6: (Color on line) Density distribution of neutrons, protons inside W-S cell for four different temperatures at ub ~ 0.050 
fm~^ and proton fraction Yp = 0.45. 
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surmounted in the regime of mean field results. 



